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We discuss a general and efficient approach for "bootstrapping" short-time correlation data in 
chaotic or complex quantum systems to obtain information about long-time dynamics and station- 
ary properties, such as the local density of states. When the short-time data is sufficient to identify 
an individual quantum system, we obtain a systematic approximation for the spectrum and wave 
functions. Otherwise, we obtain statistical properties, including wave function intensity distribu- 
tions, for an ensemble of all quantum systems sharing the given short-time correlations. The results 
■ are valid for open or closed systems, and are stable under perturbation of the short-time input data. 

f*"^ ' Numerical examples include quantum maps and two-dimensional anharmonic oscillators. 
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I. INTRODUCTION 

When a quantum system is known to have a chaotic classical limit, the simplest description of its eigenvalues, 
eigenstates, and dynamics is given by the universal predictions of random matrix theory (RMT) .1] and the closely 
i— i related random wave hypothesis Q. Recently, however, there has been increased interest in understanding the 
f""^ | deviations from RMT that are quite sizable in many systems of interest and are often due to short-time dynamical 
effects. RMT assumes implicitly that under time evolution, any initial state immediately spreads randomly over the 
entire available Hilbert space; any realistic chaotic dynamics, however, maintains at short times information about 
the initial state and only after some finite mixing time T m j x does the dynamics become truly random. 

For the purpose of describing spectral statistics, such as the distribution of level spacings, the above consideration 
may easily be taken into account by noting that RMT predictions are valid only inside energy windows of size less 
than a ballistic Thouless energy E c ~ ?i/T m i x . The situation with wave functions is not so simple, as short-time 
dynamical effects can lead to large deviations from RMT not only for energy-averaged quasimodes but also for 
individual eigenstates. 

Any short unstable periodic orbit is one obvious example of a non-random dynamical feature that is known to cause 
significant deviations from RMT statistics in the eigenstates. It has been shown that the "scarring" of individual 
wave functions by a typical periodic orbit is an O(h ) effect that persists in the semiclassical limit, as measured by 
the distribution of wave function intensities on the periodic orbit in Husimi phase space |3| ■ Using a linearization of 
the dynamics around the specific orbit, the distribution of quantum intensities on or near the orbit may be obtained 
semiclassically as a function of the classical monodromy matrix, and various moments of the distribution, such as 
the inverse participation ratio, may be expressed analytically Q. The scarring effect has been studied experimentally 
and numerically in a wide variety of systems [aL and may have consequences for the conductance through a resonant 
£h ■ tunneling diode or a ballistic quantum dot [ET3- 

The imprint of short-time dynamics on eigenstate structure is not limited to situations where the short-time dynam- 
ics is related to classical unstable periodic orbits. For example, scar-like resonances related to diffractive trajectories 
' have been observed in a 2DEG conductance experiment . The approach of using short-time behavior to understand 
eigenstate structure and statistics has been used successfully in many situations where a proper classical limit does 
not exist, such as two-body random interactions in a many-fermion system |9( or dynamics on a quantum graph or 
lattice ^3 > as wei l as m pseudo-integrable systems where the Lyapunov exponent vanishes Furthermore, short- 
time dynamical information necessarily implies deviations from RMT not only for individual wave function intensities, 
but also for spatial or phase space correlations 0. 

Our aim here is to discuss a systematic, general, and efficient framework for studying the constraints that short-time 
correlations impose on the eigenstate structure, regardless of whether such short-time correlations can be computed 
semiclassically. We allow ourselves to focus on one or an arbitrary number of initial wave packets, and the calculations 
may be performed equally well in closed or open systems, as there is no assumption of unitarity in the dynamics. We 
explicitly allow for the presence of errors in the short-time correlations, and demonstrate the stability of the results 
with respect to such errors. 

In the context of extracting stationary properties from a time-domain correlation function, we mention the important 
work that has been done by Mandelshtam and coworkers using the "filter diagonalization" method In that 

approach, one typically begins with a single wave packet, and spectral information can in principle be computed 
exactly when the correlation function is known for at least N times, where N is the Hilbert space dimension. In the 
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bootstrapping approach, we use multiple initial wave packets, and we do not assume exact finite-dimensionality of 
the Hilbert space. Thus our goal is not an exact solution of the spectral analysis problem, but rather an increasingly 
good approximation as the amount of input data increases. Of course, exact solutions are not possible in any case in 
the presence of noise or numerical instability, so in practice regularization must be performed, yielding comparable 
results for the two approaches. One advantage of the bootstrapping approach is that the linear algebra involved 
requires M x M matrices only, where M is the number of initial wave packets, allowing the problem to scale very well 
computationally for large system sizes and long times. 

The paper is organized as follows. In Sec.[nJ we define the quantities of interest in the time and energy domain, and 
obtain general expressions for the bootstrapped correlation function and spectrum. Sec. IIIII serves to define two sets 
of numerical models, which may be used to illustrate the general formalism. Next, in Sec. II VI we discuss convergence 
properties of the bootstrapping approximation, including a treatment of stability in the presence of noise. In Sec. E] 
we examine how bootstrapping may be used to compute statistical quantities of interest, including wave function 
intensity distributions and intensity correlations, using a very small amount of time-domain data as input. Finally, 
the key conclusions are briefly summarized in Sec. IVII 



II. BOOTSTRAPPING FOR CORRELATION FUNCTIONS AND SPECTRUM 



We begin by considering a set of M wave packets (pi, for simplicity taking the wave packets to be orthonormal (but 
not a complete set). In practice, the choice of <pi will be dictated by the physics of interest. For example, the 4>i 
may be chosen as position eigenstates if we are interested in position-space wave function structure, or momentum 
eigenstates for scattering behavior, or Gaussian wave packets for analyzing the effects of classical periodic orbits. In 
a many-body problem, the 4>i may usefully be taken as the non-interacting product states. 

The quantity of interest in the time domain is the correlation function 

C^-(t) = (0 j |e-^*/»|0.) j (1) 

whose diagonal elements Cu(t) constitute the autocorrelation functions for the individual wave packets. Knowledge of 
the exact correlation function for all discrete times t = mTo leads by Laplace transform to the discrete-time Green's 
function 

oo 

g 13 (e) = (iny'ToJ^ 

m=0 

« {ihy 1 / dt eWCH® = (frl , Ifr) , (3) 

Jo hi — H + te 

where the continuous-time limit of Eq. is obtained for To <C h/SE, and SE is a typical energy spread in the wave 
packets. 

We will see that it proves useful to decompose the return amplitude of Eq. Q at time t into a special "new" 
component that is returning for the first time to the subspace spanned by the M wave packets fa and the remainder 
due to propagation that has revisited this subspace at least once at time steps in between and t. In spirit, this is 
reminiscent of the T— matrix approach of Bogomolny 14], where an integral kernel is defined in terms of all classical 
trajectories that start on a given Poincarc surface of section and return to the surface of section without intersecting 
it at intermediate times. The T— matrix is defined directly in the energy domain whereas we begin our analysis in the 
time domain and transform to the energy domain later on. The decomposition used here also resembles somewhat the 
one used in the quantitative analysis of periodic orbit scarring [j| , where it is helpful to separate terms in the return 
amplitude that are associated with paths staying on the periodic orbit from terms associated with homoclinic paths 
that leave the orbit once and first return at some later time. Of course, in the case we consider here the M— dimensional 
subspace spanned by the wave packets <pi will not in general have any connection with a particular classical structure 
such as a periodic orbit or surface of section. Instead the choice of <f>i is governed either by our exact or approximate 
knowledge of the correlation function for those initial and final states or by an interest in extracting wave function 
structure information in a specific basis or phase space region. Furthermore, no semiclassical approximation is implicit 
in the method we develop here, although we will see below that the approach can be extended to situations where 
the correlation function information used as input is only approximate, as would be the case for example when a 
semiclassical propagator is used. 

Formally, the new recurrences Bij (m) at time t = jtiTq may be defined as 

%(m) = Me-** 1 '''/* ( (1 _ py-iHTo/ny- 1 ^ (4) 
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where 

M 



P = J2\^)(M (5) 



fe=l 

is the projection onto the subspace of interest. More explicitly, these new recurrences may be computed from the full 
Cij amplitudes as 

o / x JC«(To) rn = l 

^ (m, -\C l3 (mT )~j:™=ii:lLiB t k(p)C kj ((m-p)T Q ) m>2 ' W 

The full correlation function is then given by a convolution, 

m A/ 

a, (mT ) = ^ 2 B lfe (p)C*j ((m - p)T ) (7) 
p=i fe=i 

or, in matrix notation, 

m 

C(mTo) = 53B(p)C((m-p)T ) 

P =i 

m— 1 m — 2 m— p— 1 

= B(m)+53B(p)B(m-p)+53 £ B(p)B(|/)B(m - p - p') + • ■ ■ , (8) 

p—i p— l p'_ i 

where C(0) is always the identity matrix. The matrix B(m) records amplitude that at the m— th step returns for the 
first time to the subspace spanned by the <pi, while terms in Eq. |(HJ involving a product of n B-matrices correspond 
to processes where amplitude leaves and returns n times to the same subspace over m steps. In the energy domain, 

Gail?) = (tfi)- 1 !^ ^* - , (9) 

where 

oo 

B(E) = 53 e mToE/h B(m) . (10) 

m=l 

We are however interested in the information that can be extracted from knowledge of the correlation function at 
a finite set of times t only, say t = tuTq for m = 1 ■ • ■ L, possibly in the presence of noise. If we assume C(mTo) is 
known only for times t < T max = LT$, i.e., for 1 < m < L, then we may compute the new recurrences B(mTo) only 
for 1 < m < L using Eq. ||BJ|. It is convenient to define 



„ , v f B(m)e- mT °/ T 1 < m < L 
' v ; [0 otherwise 



(11) 



The cutoff time r, which can be much larger than the bootstrap time T max = LTq, serves as a smoothing scale in 
the energy domain, and its significance will be discussed in Sec. IIVI below. Given just the matrices Bi. T (?n), we may 
compute a "bootstrapped" approximation to the full correlation function at all times: 

m 

C l , t (toT ) = 53b l , t (p)C l , t ((to-p)T ) 

p=i 

m— 1 m—2m—p—l 

= B L , r (rn)+ 53B L:T (p)B L , r ( ? 7i-p)+ 53 53 B L:T (p)B L , r (p')B L , T (m - p - p') + ■ ■ ■ , (12) 
p— l p—i p'—i 

having the property Cl jT (?tiTo) = C(mTo)e~" iT °/ T for m < L. In the energy domain, 'Bl^(E) may be defined as a 
Laplace transform of Bi lT (m), in complete analogy with Eq. Hl()(l above. 
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Finally, one often encounters a "noisy" situation where even the short-time dynamics is only approximately known. 
For example, we may be interested in building up the full dynamics using only semiclassical expressions for the 
propagator at short times. We then have knowledge of 

C e (mT ) = C(toTo) + eD(mTo) (13) 

for 1 < m < L 1 where D(mT ) are quasi-random, uncorrelated error matrices and e characterizes the relative size 
of the error. Given this input data we may calculate approximate "new" recurrences Bi ir c (mTo) by extending the 
exact formula of Eq. (jSJ, 

r C e (T ) e - T «/- m=l 
Bi, T>£ (m) = I C e {mT )e- mT °/ T - Y^=i B L , r , e (p)C e ((m - p)T )e-( m -P^>/ T 2 < m < L . (14) 
[ otherwise 

The "bootstrapped" long-time evolution Cz, jTi£ is given by iterating these approximately known short-time "new" 
recurrences analogously to Eq. l(T2*|l . 

rn 

C L ,rA mT o) = ^ B L . T . e (p)C L . T . e ((m - p)T ) 

p=l 

m—1 m—2m—p—l 

= B L ^ e {m)+J2 B L.r.Ap)' B L.r.A m -p)+J2 B L,r,e(P) B L,r,e (p') B L,r,e (m - P - p*) + • <15) 

p=l p=l p' = l 

Again, by construction the bootstrapping procedure simply reproduces the noisy input data for times t below the 
bootstrap time T max , i.e., Cl.t.h = C e e -mT °/ T for m < L. However, bootstrapping allows us also to learn something 
about longer times t > T max using the short-time correlation function. 



III. NUMERICAL MODELS 



A. Quantum Maps 

Classical and quantum chaotic maps in one dimension are often used as the simplest examples for illustrating general 
chaotic behavior, and share many scaling and other physical properties of two-dimensional Hamiltonian dynamics |15| . 
Discrete-time maps may be thought of as arising from a continuous-time Hamiltonian dynamics either via a Poincare 
surface of section or by stroboscopically viewing motion in a periodically driven Hamiltonian. In the latter case, we 
may consider 

^ oo 

H(q,p,t) = -=—T{p) + V(q) V 5(t - jT kick ) , (16) 

J kick . ' 

which yields 

Pj+i = Pj-V'(qj) 

qj+i = qj + T'(p j+1 ) , (17) 

when the position qj and momentum pj are recorded just before kick j. The corresponding quantum evolution over 
one time step is given by 

U = e -iT(P)/n e -iv( 9 )/n (18) 

As a specific example, we may take T(p) — \w p p 2 -I- K p (s\a.p — i sin 2p) and V(q) = — \w q q 2 — K q (sin q~\ sin 2cj) , 
for a toral phase space (q,p) G [— vrO, tt) x [—it, it). With integer values of w q and w p , this is a perturbed cat map 16] : 

Pj+i = pj + w q qj + K q (cosqj — cos2qj) mod 2ir 

Qj+i = Qj + w pPj + Kp (cos pj — cos 2pj) mod 2ir , (19) 

where nonzero values of K q , K p are needed to break the symmetries and ensure nonlinearity of the dynamics. One 
easily checks that the dynamics is completely chaotic for sufficiently small K q ^ p . 

For this compact classical phase space, the quantum evolution of Eq. i|18fl acts on a Hilbert space of dimension 
N = 2tt/H : the mean energy level spacing is A = 2Trh/NTy iic ] i = h 2 /T^ck, and the Heisenberg time at which levels are 
resolved is Tjj = NT^^. Since the map dynamics is already discretized, it is natural to use the period as the 
time step Tq in the bootstrapping calculation. Without loss of generality, we may choose units where T$ = Tkick = 1- 
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B. Two-Dimensional Wells 



As our model of a non-integrable system with a time-independent Hamiltonian, we use the Barbanis Hamilto- 
nian |l7j . which describes a two-dimensional anharmonic oscillator: 

H(x',y',p' x ,p' y ) = + + I mw >' 2 + + Xx'y' 2 . (20) 

After a canonical transformation and an overall rescaling of the energy, the Barbanis Hamiltonian may be re-written 
as 

H(x, y ,p x ,p y ) ^^ + P l+ 1 -^ + ^ + |xy 2 , (21) 

where a is a dimensionless parameter characterizing the shape of the well. In these dimensionless coordinates, the 
metastable well has its minimum at x = y = and extends from x = — 1 to x = 1 along the y — symmetry axis; 
the barrier height is E max = 1/2. Upon quantization, one additional parameter besides a is introduced, namely an 
effective Ti or equivalently the number of quantum levels below -E max . 

As compared with the simple quantum map model presented above, analysis of bootstrapping in the Barbanis 
system requires consideration of the following three circumstances, which are typical of many Hamiltonian systems: 
(i) time is not naturally discrete and thus an explicit choice is needed for the discretization time To, (ii) resonances in 
the metastable well have finite intrinsic width, introducing a new long-time scale, and (iii) classical dynamics in the 
well is mixed, with the phase space at energies of interest shared by regular islands and a chaotic sea. The implications 
of these three circumstances will be discussed below, when numerical results for bootstrapping in the Barbanis system 
are presented in Sec. IIVBI 



IV. CONVERGENCE PROPERTIES AND SENSITIVITY TO ERROR 



In this section, we examine how bootstrapping may be used when the given information about the short-time 
correlation function is sufficient to compute (approximately) the long-time dynamics and spectrum. An alternative 
situation, where the given information only restricts us to an ensemble of possible long-time behaviors, and the 
objective is to obtain statistical properties of the long-time dynamics or spectrum, is discussed in Sec. 



A. Noise-Free Bootstrapping 

We want to estimate the error made in using short-time information up to the bootstrap time T max = LTq to 
estimate long-time dynamics in a chaotic system at times t ^> T max . Let us first assume negligible noise by setting 
e = in Eq. (|13l) . Clearly the error is then associated with amplitude that starts in the subspace spanned by the M 
wave packets <pi and is never captured by the short-time correlation function because it does not return to the original 
subspace at any time during the first L steps of evolution. In terms of the B— matrices discussed in the previous 
section, the total probability that does not return in time r max = LTq is 



1 L 

P(T max ) = 1 - — ^ B(m)tB(m) . (22) 



M 

m—l 

Mathematically, the probability P(T max ) is clearly related to the probability of remaining for at least time T max in 
a system with M maximally coupled open decay channels. When the dynamics is chaotic, this probability can be 
represented analytically as an integral in the context of random matrix theory [l^ ; for our purposes it is sufficient to 
note that 

( p -MT m .^/T H J- ^TtjIJM 

r l J mai — \ 1 (T IT 1 \M+1 rp „ rp > l ZJ j 

I M+T y 1 H I 1 ^maxj Jmax^JiJ 

where Tjj is the Heisenberg time, and the power-law long-time limit also serves as an upper bound for P(T max ). 
Clearly, we require T max > Th/M in order to recapture most of the initial amplitude, so that the lost probability is 
small. We emphasize that this estimate, based on random matrix theory, may be used to obtain the correct scaling 
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FIG. 1: The lost probability P(T max ) associated with eliminating new recurrences arriving after time T max is plotted as a 
function of T max for the quantum map of Eq. I|19|l . with w q = w p — 1 and parameters K q , K v uniformly distributed between 
— 1/2 and 1/2. Squares, circles, and triangles represent data when the number of wave packets M is M = 1, 2, or 4, respectively. 
Open symbols are associated with system size TV = 64 and closed symbols with system size N = 128. The dashed curve is the 
small- T max (classical) limit in Eq. i'l'M . while the three solid lines represent the large- T max power-law falloff for M — 1, 2, and 
4. All quantities shown in this and subsequent figures are dimensionless. 



behavior of the lost probability P(T max ) with bootstrap time T max , even when the prefactor in Eq. Q23[) is invalid due 
to nonrandom short-time dynamical effects. 

The behavior of the lost probability P(T max ) for small and large T max is illustrated in Fig. ^ Here an average 
over quantum maps given by Eq. (|19|1 has been performed, with w q = w p = 1 and nonlinearity parameters K q , K p 
randomly distributed between — 1/2 and +1/2. We note the expected exponential behavior for small T max , with the 
classical decay rate M/Tjj, as well as the rapid power-law decay of the lost probability for T max > Tjy, especially in 
the case of multiple wave packets M > 1. 

We are interested in the error induced at long times t 3> Tmax by omitting new recurrences not captured in the 
short-time correlation function. At this point, we have not introduced any smoothing of the input data, i.e., r = oo 
in Eq. (|ll(l . The typical returning amplitude at time t has completed 0(Mt/Tu) cycles of leaving and returning to 
the subspace spanned by the M wave packets 4n, i.e., in Eq. (jSJ) the dominant terms are ones involving a product of 
0(Mt/Tji) B— matrices. In each cycle, probability given by Eq. (|23|l is lost, with the errors accumulating coherently, 
so that the relative error in the matrix elements at time t 3> T max = LTq > Th/M is given by 

M-C^ /MA 2 

m) - \mw l?v P(Tmax) ' (24) 

where ||C(i)|| 2 = Tr C^C^) = £V |Cy(<)| 2 - The quadratic growth in the long-time error is clearly seen in Fig. [21 
for several choices of the bootstrapping parameters. 

To find the time scale Tbrcak beyond which the bootstrapping procedure breaks down, we assume Th/M < T max < 
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FIG. 2: The ratio of the relative error in the propagator, E{i), to the probability lost on average during each cycle of recurrences, 
-P(Tmax), is plotted for times t > T max , M = 2, several values of T max , and two system sizes. All system parameters are the 
same as in the previous figure. Squares, circles, and triangles correspond to MT^x/Th = 5, 10, and 20, respectively, while 
open and closed symbols distinguish system size N — 64 from system size N = 128. The solid line indicates the quadratic 
growth of the error consistent with Eq. I|240 ■ 



Th I \/M. Then, setting the right hand side of Eq. (|24J) to unity and using Eq. Q23|l , we obtain 

J break ~ 1 max , ™ /rr ■ 

We see that including only a minimal number of new recurrences by setting T max ~ Th/M leads to breakdown 
of the bootstrapping approximation soon thereafter (Tbreak ~ T max ), but including additional recurrences leads to 
exponential growth in the accuracy of the bootstrapping approximation and consequently to exponential increase in 
the breakdown time. Of course, this exponential growth ceases at very large values of T max , when the error becomes 
dominated by a small fraction of eigenstates that have unusually little overlap with the wave packets fa. Then P(r m ax) 
follows the power-law behavior of Eq. (|23|l , and the growth in Tb rea k accordingly crosses over to a power-law behavior 
with T max : 



T brca k~T ff ^M^(^) ,._>,„ 



M \ T, 



<-H 



for T max > Th- We note that the growth of the break time Tb roa k with increasing bootstrap time T max remains faster 
than linear except in the single- wave packet case M = 1. This super- linear growth is illustrated in Fig. [3] for the case 
of two wave packets (M = 2), where the break time Tbreak has been quantified as the time scale where the relative 
error of Eq. I|24(l reaches unity. 
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FIG. 3: The break time Tb rea k of the bootstrapping approximation, defined by -B(Tb re ak) = 1, is plotted for M = 2 and several 
choices of Tmax. Open and closed squares represent system sizes N = 64 and N = 128, respectively. The ensemble of systems 
is the same as in the previous two figures. The straight line is the theoretical prediction Tbreak ~ TmaL +1 ^ 2 ~ Tmix of Eq. 1261 . 

B. Results in the Energy Domain 

Starting with known short-time information about the correlation function, the bootstrapped long-time dynamics 
may be Laplace or Fourier transformed into the energy domain to obtain good approximations to the Green's function, 
spectrum, or local density of states. Alternatively, the short-time "new" recurrences may be transformed directly into 
the energy domain to obtain spectral information, as indicated by Eqs. @ and (|1U|) . To avoid unphysical oscillations 
in the spectrum on energy scales below 7i/Tbrcak (associated with the breakdown of the bootstrapping approximation 
at long times), we impose an explicit smooth cutoff on the the short-time dynamics, in accordance with Eq. I|ll[) . Loss 
of information is minimized by choosing the cutoff time r of the order of Tb roa kj which is equivalent to Lorentzian 
smoothing of the spectrum on the scale ft/Tb rea k. 

The numerical data presented in Fig. 0] is obtained for the Barbanis potential, with parameters a = 1.1 and 
h = 0.0198 in Eq. (|21|l . corresponding to slightly over 300 quantum resonances in the metastable well. Six initial 
Gaussian wave packets are used in the calculation, all centered at x = and having an average momentum of 
magnitude \p\ — 0.96, corresponding to an energy E — 0.46. Thus, we are viewing dynamics slightly below the top 
of the barrier, E max = 0.5. The classical dynamics in the energy range considered is approximately 57% chaotic, as 
measured using a Poincare surface of section at x = 0. All six initial wave packets are centered in the chaotic region 
of classical phase space. 

We see from the middle solid curve in Fig. 0] that most peaks in the spectrum can readily be resolved using 
bootstrapping, taking correlation information up through the Heisenberg time Tjj as our only input. For bootstrap 
time T m ax = 2Tfj (top solid curve), the spectral peak heights already stand out by four orders of magnitude above 
the background. The root mean squared error in the bootstrap-predicted peak locations drops from 0.034A when 
?max = Tjj to 0.0064A for r m ax = 2Tff, where A is the mean level spacing. We may contrast this with the result, 
indicated by the dashed curve, of merely transforming and smoothing the same correlation information, up through 



9 



10 16 




-j qu l ^ i ^ l , i , : i : i , 

0.452 0.454 0.456 0.458 0.46 

E 

FIG. 4: The local density of states summed over M — 6 wave packets located at energy E — 0.46 in the Barbanis potential 
of Eq. ED, S(E) = J2 =i Re (»/"■) Gjj(E), is computed using the bootstrapping approximation in accordance with Eq. JJjJ. 
The three solid curves from bottom to top correspond to different bootstrap times used in the bootstrapping calculation: 
Tmax = Th/2, T m ax = Th , and T max = 2Th ■ For comparison, the result of Laplace transforming the correlation function 
through time T max = 2Th, without bootstrapping, is shown as a dashed curve. Each spectrum has been scaled by an arbitrary 
constant to allow for easy comparison on a single plot. In each case a smoothing time scale r ~ Tb re ak has been chosen to 
remove unphysical oscillations in the spectrum in the energy range shown. The dotted vertical lines indicate the locations of 
the exact resonance peaks, obtained by taking T max — > oo. 



Tmax = 2Tff, but without the benefit of bootstrapping. Here the resolution is very poor, and we are far from being 
able to detect, for instance, the two doublets near E = 0.4515 and E = 0.4535. 

In contrast with the map model studied in Section IIV Al in the Hamiltonian system investigated here we must 
discretize time explicitly by introducing a new time scale To. The results of the calculation, however, are unaffected 
by the choice of To, as long as Tq < h/SE, where SE is the energy uncertainty in the wave packets <pi. Equivalcntly, 
the time step T must be chosen short enough so that the self-overlaps Cu(T ) are large due to free-flight dynamics. 

A second key difference with the map model is that quantum motion in the Barbanis potential is described by 
resonances rather than bound states. Indeed, by comparing the upper two curves in Fig. 0] we see that in the 
Tmax = Tjj bootstrapped spectrum, the widths of several peaks (e.g., the rightmost one near E — 0.4612) are already 
dominated by the intrinsic resonance widths rather than by any error associated with the time cutoff. In general, the 
efficiency of the bootstrapping approach increases as one considers systems that are more open, since it is sufficient 
to choose a bootstrap time T max that will generate accurate dynamics to time Tb rea k ~ ^decay, where Tdocay is the 
intrinsic lifetime of the resonances, possibly shorter than Th- 

Finally, a third major difference between perturbed cat maps and the Barbanis potential is the presence of regular 
as well as chaotic states in the Barbanis spectrum. By choosing the test wave packets (pi appropriately, one may 
optimally resolve states in the phase space region that are of greatest interest in a given application. For example, 
in ordinary scar theory, one may begin with a wave packet centered on a specific periodic orbit Q, with the aim 
of obtaining optimal information on the structure of wave functions with high intensity on that orbit and their 
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associated eigenvalues; the price to be paid is the suppression of the "anti-scarred" eigenstates that have anomalously 
low intensity on the same orbit. Here, we have randomly placed the six test wave packets in the chaotic portion of 
phase space, improving our ability to resolve the chaotic states, but necessitating the use of longer bootstrap times 
T max to identify spectral peaks associated with regular wave functions, such as the very narrow resonance peak near 
E= 0.4536. 

In this context we note also that, in the case of scar theory, little or no benefit is gained by following several wave 
packets launched along the same weakly unstable periodic orbit, since they all exhibit very similar time evolution, 
and share nearly identical local density of states |19j . From a bootstrapping perspective, we may consider two wave 
packets nearly related by time evolution, e.g., e~ %Hr i n \(f>\) . Then the probability 1 — P(T max ) of Eq. 122|) 

for returning to the subspace spanned by 4>\ and <f>2 by time T max is nearly the same as the probability of returning 
to 0i itself, assuming T max 3> r. The rapid decrease in the "lost probability" P(T max ) with increasing number of 
wave packets M, as indicated by Eq. depends entirely on the wave packets behaving in an uncorrelated manner. 
Thus, the bootstrapping procedure for multiple wave packets is most effective when the wave packets are chosen from 
different regions of phase space to avoid obvious correlations. 



C. Influence of Noise 



Noise in the input signal may be an important factor in specific applications of the bootstrapping procedure, for 
example where a semiclassical or other approximation is used to calculate the short-time correlation function. We 
return to the quantum map model discussed in Sec. IIV Al and introduce random noise into the short-time correlation 
function, as indicated in Eq. 1|13|) . The random error matrix elements Dij(m) in Eq. I|13[l are independent Gaussian 

random variables of zero mean and variance 1/N, where N is the Hilbert space dimension, so that \Dij{m)\ 2 = 
\C'ij(m)\ 2 at long times m. Then the dimensionless parameter e characterizes the relative size of the noise. A 
spectrum may be produced by bootstrapping the noisy short-time data. The results of such a calculation are presented 
in Fig. |SJ We see that the spectral reconstruction is quite robust for small e, and breaks down at around e = 0.2 or 
0.3, independent of N and M. 

To study more carefully the breakdown in the accuracy of the bootstrapped spectrum and its dependence on 
parameters N and M, we need to define a quantitative measure of the error in the bootstrapped spectrum. Consider 
a local density of states S(E) reconstructed from the exact correlation function known for t < T max and a local density 
of states S e (E) reconstructed from the same input but with added noise characterized by e as in Eq. ljT3|) . We may 
define the dimensionless error ratio 

Z[e) - JdE (In S(E)r ' { ' 

which measures the error induced in the reconstructed spectrum by noise of size e in the input. We note that it is 
appropriate to focus on the logarithm of the reconstructed spectrum, because the spectrum itself is dominated by 
sharp peaks, as seen in Fig. [S] The quantity Z(e) is shown in Fig. for the same parameters as were used earlier 
in Fig. Not surprisingly, we observe growth in the spectral error with increasing noise e, but, more importantly, 
this error is almost independent of system size iV and number of wave packets M (at e = 0.3, Z(e) varies at most 
by 30% as N changes by a factor of 4 and M by a factor of 3). The same results have been observed for other 
system parameters. This robustness implies that input with noise of a small but finite size e may be used in the 
semiclassical limit N — > oo (equivalently, h — > 0). Obviously, an even more favorable situation exists when the noise 
level e decreases with increasing N. An important example of the latter situation exists when the "noise" results from 
using a semiclassical (stationary phase) approximation for the short-time dynamics t < Tmax- Then e ~ h, which 
decreases with N. Therefore short-time evolution computed within a semiclassical approximation may be used with 
confidence to obtain stationary properties of the exact quantum system. 



V. BOOTSTRAPPED WAVE FUNCTION STATISTICS 



We now consider the situation where known short-time dynamical information is insufficient for resolving individual 
eigenstates, and the focus therefore shifts to predictions of a statistical nature. In other words, we consider an 
ensemble of systems that all share (perhaps approximately) a given short-time dynamics, and ask what information 
can be extracted about the distribution of wave functions in systems drawn from this ensemble. 
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FIG. 5: The local density of states summed over M randomly placed wave packets, S(E) = ^2 = i R e (V 71 ") Gjj(E), is computed 
for the quantum map of Eq. 1191 . with w q = w p — 1, K q = 0.2, and K p — —0.3, using the bootstrapping approximation after 
noise has been introduced into the short-time input data. From top to bottom, the three sets of curves correspond to (a) 
N = 128, M = 2, (b) N = 128, M = 6, and (c) N = 512, M = 2, where N is the system size or Hilbert space dimension. 
Within each set, the top (solid) curve is the reconstructed spectrum in the absence of noise, and the three dashed/dotted 
curves, from upper to lower, indicate reconstructed spectra for the same system with the dimensionless noise parameter set to 
e = 0.1, e = 0.2, and e = 0.3. Each spectrum has been scaled by an arbitrary constant to allow for easy comparison on a single 
plot. In all cases, the spectrum is reconstructed from the correlation function for t < T max = 2>Th/M. 



A. Inverse Participation Ratio Calculations 

The simplest quantitative measure of wave function structure is the inverse participation ratio (IPR) or second mo- 
ment of the wave function intensities [2(|: IPR* = N^2f =1 | (0.;| 1 I') | 4 , where ^ is an eigenstate, TV is the Hilbert space 
dimension, the sum is over a complete basis and we impose the usual normalization condition J^., |(</>i| , J , |)| 2 = 1. 
The IPR measures the degree of wave function localization, ranging from 1 for a delocalized wave function having 
uniform overlaps with all basis states to N for a completely localized state v]>. RMT predicts IPR = 2 in the ab- 
sence of time-reversal or other symmetry. Similarly, for each basis state <fi we may define a local IPR (LIPR) as 
LIPR = iV£}=i | | 4 , where the sum extends over eigenstates Q ; the LIPR measures the degree of localization 
associated with a specific basis element <fi and is proportional to the average long-time return probability |(</>|</)(i))| 2 
as t — > oo. 

Extending arguments developed originally for periodic orbit scars 0, , we may interpret long-time dynamics in a 
chaotic system as a convolution of known short-time recurrences with quasi-random long-time recurrences, 

T 

(M(t)) « {4>\<t>{r))r^{t-T), (28) 

T=-T 

where for simplicity we have assumed discrete-time dynamics, the sum over r extends to some scale T that includes 
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FIG. 6: The error Z(e) in the reconstructed spectrum for the quantum map of Eq. 11911 is computed as a function of the noise 
level e in the input correlation function. All parameters are identical to those in Fig. The solid curve indicates system size 
N — 128 with M — 2 wave packets, the dashed curve is for N = 512 with M — 2, and the dotted curve is for N = 128 with 
M = 6. In all cases, the noise-free spectrum S(E) and the noisy spectrum S € (E) are both reconstructed from the correlation 
function for t < T max = 3T H /M. 



as much as possible of the non-random dynamics of interest but is still short compared with the Heisenberg time 
Th, and r$(t') are Gaussian random independent variables, associated with nonlinear long-time recurrences. For the 
LIPR, we obtain 



LIPR, 



T 

! E 

r=-T 



(29) 



(30) 



in the notation of Sec. [H] where the overall prefactor of 2 is the RMT result in the absence of time-reversal sym- 
metry. The autocorrelation function or return amplitude C^(r) may be computed from the B-matrices using the 
bootstrapping formulas of Eqs. JBJ and JJJ, or we may explicitly write 



LIPR, 



OO 

r=l 



T m ax 

E 



-B^(ti)<5(ti - t) 



E b 00'( t 2)-B0'0( t i)^(ti + r 2 - t) 

Tl,T2 = l 0' 



(31) 



Here the upper limit T in the sum over r may safely be taken to infinity, as long as the bootstrap time Tmax Th/M, 
so that most of the probability is lost by the Heisenberg time, and times r ~ Th do not contribute significantly to 
the sum. We note that the second- and higher-order bootstrapping terms implicitly include revivals at times longer 



13 



than T max , although only the correlation function up to T max is used as input to the calculation. The bootstrapping 
formula makes optimal use of the available short-time information, and good agreement may be obtained even for 
fairly short bootstrap times T max . 

As a specific example, we consider a quantum map defined by Eq. I|17|) . with kinetic term T(p) = ip 2 and kicked 
potential 



V(q) 



1 



(q-qo? 



-&(q~qo) + ^-^e(qo-q) 

qo i - qo 



(32) 



where as before q and p both range from — 7r to +7r, and 0(x) is the usual step function defined by Q(x) = 1 for 
x > and 6(2;) =0 otehrwise. The dynamics is fully chaotic and has no period-one classical orbits, but does have 
a diffractive orbit at q = qo, p — 0, associated with a cusp in the kick potential. The bootstrapping calculation is 
performed for a single wave packet centered on this diffractive orbit. In Fig.0 we calculate the LIPR for this wave 
packet, as a function of parameter vq, exactly and in the bootstrapping approximation. We see immediately that 
RMT (equivalent to the bootstrapping prediction with r max = 0) severely underestimates the degree of wave function 
localization when vq < 1 and the diffractive orbit is consequently strong. Bootstrapping the one-step recurrence only 
(Tmax = 1) greatly overestimates the effect, but the T max = 2 calculation, which incorporates information about 
one-step and two-step new recurrences, already gives a good approximation to the exact answer over the entire range 
of vo. We note that the bootstrapping has been performed here using one- and two-step time correlation data for a 
single wave packet; obviously the results can only improve if multiple wave packets are used simultaneously with the 
same bootstrap time T max . 
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FIG. 7: The local inverse participation ratio (LIPR) at the location of the cusp is calculated for the diffractive potential of 
Eq. 132L for system size N = 64, one wave packet centered on the cusp at qo = — 0.2tt, and several values of the kick parameter 
Do- The exact data is averaged over boundary conditions for each value of vo. The bootstrapping prediction, using Eq. (IMOli 
for one wave packet centered on the cusp, is shown for T max = 1 and 2, in units of the kick period. The RMT prediction, 
equivalent to bootstrapping with T max = 0, provides a baseline for comparison. 

We now fix vq — 0.59 and repeat the above bootstrapping calculation for single wave packets centered at various 
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locations in phase space. In each case, we find the exact LIPR by diagonalizing the evolution matrix. We also predict 
the LIPR using the bootstrapping approximation with T max = 3, i.e., the recurrences for three time steps are computed 
exactly, bootstrapped to obtain long-time behavior, and then used to estimate the local inverse participation ratio in 
accordance with Eq. i|3U|) . The results are shown in Fig. [S] Here, the bright spot slightly to the left of center is the 
localization peak associated with a diffractive orbit at q = qo = — 0.2ir, p = 0. We observe that the bootstrapping 
procedure allows not only this peak but most significant features of the localization landscape to be well resolved by 

max u • 
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FIG. 8: The local inverse participation ratio (LIPR) is plotted as a function of position and momentum for the potential of 
Eq. i'-')'2t . with system size N = 64, cusp location go = — 0.27T, and cusp strength «o = 0.59. The exact LIPR landscape is shown 
in the upper panel, while the lower panel represents the prediction of a bootstrapping procedure with T max = 3. The color 
scale ranges from LIPR=2 (black) to LIPR=5 (white). 



B. Wave Function Intensity Distribution 

A prescription similar to the above may be used to compute higher moments of the intensity distribution beyond 
the standard inverse participation ratio; instead, we turn our attention to the intensity distribution itself. Knowledge 
of such a distribution is essential, for example, to the understanding of resonance width or decay rate statistics in a 
weakly open system. In the context of scarring, it has been shown that the probability distribution of wave function 
intensities may be obtained by combining a smooth spectral envelope constructed from the short time dynamics 
with Gaussian random fluctuations on fine energy scales [|lj . More generally, whenever a separation of scales exists 
between non-random short-time dynamics and quasi-random long-time behavior, we may write the local density of 
states (strength function) for wave packet o as [jj 

Re (i/n)G^(E) = £ 6(E - E n )\{^ n )\ 2 w £ S(E - E n )S s ^° oth (E)\R n \ 2 , (33) 

n n 

where 

T 

Sl mooth (E) = {^W))e lEt,h (34) 

T=-T 

is a Fourier transform of the short-time signal and R n are independent Gaussian random variables with variance 
l/N (real or complex depending on the presence or absence of time reversal symmetry, respectively). The above 
expressions assume no symmetry, with the possible exception of time reversal, and must be appropriately modified in 
the presence of such symmetries, including parity invariance 4J. The multiplication in Eq. (|33|l of a known short-time 
signal by a long-time signal assumed to be quasi-random is the energy-domain counterpart of the convolution formula 
appearing in Eq. J2HJ- 

In the bootstrapping context, we may obtain the short-time envelope using Eqs. 10 and I|1U|) . where exact new 
recurrences B(m) are replaced by Bi. T (m) as defined by Eq. Ijllfl for some choice of T max = LTq and a smoothing 
time scale r. This is the same procedure we used to construct approximate bootstrapped spectra in Sec. IIV Bl except 
that there the bootstrap time T max was chosen sufficiently long to resolve individual states, T max > Th/M, while here 
we may take T max to be only a small multiple of the one-step time To. 

Once a short-time local density of states envelope S^ nooth (E) is known, we may directly construct the probability 
distribution of wave function intensities I = \(^)\^ n )\ 2 . We need only to multiply the envelope heights S^ nooth (E), 
with uniformly distributed energies E, by random factors |i?| 2 where R is Gaussian-distributed (and |i?| 2 is therefore 
exponentially distributed for complex R): 

P W=r W dE d(\R\ 2 )e~W 2 6(l-S;™° th (E)\R\ 2 ) . (35) 

Typical examples of the resulting intensity distribution are shown in Fig. 03 Here we use the same system and 
wave packet location as in Fig. [7| but fix kick parameter vq at the value 0.2. The short-time envelope S 1 ™ 100 * 11 is 
constructed either using only one-step new recurrences (bootstrap time T max = Tq — 1) or using one- and two-step 
new recurrences (bootstrap time T max = 2Tq = 2). The T max — 1 short-time envelope already results in a predicted 
intensity distribution that is a great improvement over the RMT prediction, correctly predicting an excess of very 
large and very small wave function intensities at the cusp. The T max = 2 envelope predicts an intensity distribution 
that is in even better agreement with actual data. 

C. Wave Function Correlations 

The bootstrapping approach lends itself naturally to the consideration of observables beyond the statistics of 
individual wave function intensities / = | \ 2 ■ As a simple example, we may consider the covariance COY < j >1 .^ )2 = 
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FIG. 9: The distribution of wave function intensities / = |(<^| 1 I'„}| 2 for the kicked map of Eq. 13211 with system size N = 64, 
cusp location qo = —0.2n, and kick parameter vq = 0.2 is shown, where is a Gaussian wave packet centered on the cusp at 
q = qo, p = 0. The solid curve shows the exact data, obtained by diagonalizing time evolution matrices, and averaging over 
systems with different boundary conditions. The bootstrapped predictions are computed using Eq. 1351 . where the smooth 
envelope is obtained from bootstrap time T max = 1 or T ma x = 2, for a single wave packet M = 1. The RMT prediction, 
equivalent to bootstrapping with T max = 0, provides a baseline for comparison. 



Nj2j=i l(0i|^j)| 2 | v&l^j)! 2 - 1) where <f>\ and (f>2 are two wave packets and the sum is once again over the eigenstates. 
Obviously the covariance is a generalization to two wave packets of the local inverse participation ratio discussed 
earlier: LIPR^ = COV^.^ + 1. The covariance or correlation between wave function intensities at two points is clearly 
important, for example, for understanding the statistics of conductance peak heights in a weakly open quantum dot 
with two leads pj; it is also relevant for analogous reaction rate calculations or for the computation of interaction 
matrix elements. 

Letting \<f>) — ^gd^i) + e ie \4>2)), using Eq. I|29fl for LIPR^, and averaging over the relative phase 0, we obtain 

T 

cov,^ w Yl [K^\Mr))\ 2 + <MMmMt)\<h)]-i 

T = -T 

T T 

= |C 0102 (O)| 2 + 2^|Q 102 (T)| 2 + 2Re^Q 101 (r)C; 202 (r) (36) 

T=l T—l 

Two types of terms are present in Eq. (|36|) : ones associated with the short-time probability for evolving from state 
<j>\ to state 4>2 or vice versa, and ones associated with a correlation between the individual short-time autocorrelation 
functions for <j>i and <p2- Once again, the correlation functions C$#>(t) may be computed using the bootstrapping 
formula given by Eq. Q or Eq. |JS}, where the "new" recurrences B(m) are known up to the bootstrap time T max = 
LTq, as in Eq. I|ll|) . As in the LIPR calculation, the upper limit T of the sum in Eq. Q36fl may be taken to infinity, as 
long as Tmax < Th/M. For the covariance calculation, it is most convenient to perform the bootstrapping with just 
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M = 2 initial wave packets: (pi and (f>2- 

As an example, we consider another quantum map, defined by Eq. (|17f) with kinetic term 



T{P) = 2^~Po) 2 + b sin2(p-p ) - -sin4(p-p ) 



(37) 



and periodic kick 



V(q) 




-ir < q < -7r/2 
-7T/2 < q <q 
q <q< tt/2 

7T/2 < ^ < 7T 



(38) 



This potential has a cusp-like maximum of height vq at q — —tt/2 and another of height 3vq/2 at q = tt/2, resulting 
in the possibility of diffractive periodic motion between (q = —n/2,p = po) and (q = tt/2,p = po). We compute 
the covariance between wave function intensities |(0i| x I , )| 2 and |(</>2| V I')| 2 , where 0i and 02 are Gaussian wave packets 
centered at the two points in phase space. The results are presented in Fig. 1101 as a function of the cusp height 
parameter vq. Once again, the bootstrapping prediction is shown for bootstrap time T max = 1 or 2, in units where 
the kick period T^ i( ± is set to unity. The RMT prediction, corresponding to bootstrap time T max = 0, is shown for 
comparison. Just as in the LIPR and intensity distribution calculations, rapid convergence is observed with increasing 
Tmax, and almost all relevant information is already obtained by bootstrapping the one-step and two-step dynamics. 



Short-time dynamical information, either of classical origin or otherwise, inevitably leaves its imprint on the long- 
time behavior and stationary properties of a quantum system. The bootstrapping approach allows this information 
to be processed systematically, for one or an arbitrary number M of initial wave packets. Because multiple iterations 
of the known short-time dynamics are included, the resulting spectral accuracy can be much greater than what one 
would obtain, for example, by a simple Fourier transform of a short-time signal. At the same time, the procedure 
is extremely efficient, requiring at each energy linear algebra operations involving only M by M matrices, and 
independent of the total size N of the Hilbert space. There is no assumption of unitarity in the dynamics, and the 
procedure works equally well for closed or open systems. Robustness to errors in the short-time signal implies, for 
example, that reliable calculations can be performed when the short-time correlations are computed in a small-?i or 
other approximation relevant to a given problem. 

The bootstrap time T max can be varied to extract maximum information from the least amount of input data. At 
small values of T max , the approach can be viewed as a generalization of standard periodic orbit scar theory, leading 
to statistical prediction beyond RMT for local density of states and wave function statistics. Reliable quantitative 
predictions can be obtained for inverse participation ratios, wave function intensity distributions, and wave function 
correlations, even when the short-time dynamics is of non-classical origin. Increasing either T max or M allows for a 
systematic inclusion of additional correlations. Once the product MT max becomes comparable to the Heisenberg time 
Th, it becomes possible to go beyond statistical predictions to resolve individual eigenstates and energy levels, with 
an accuracy scaling exponentially with MT max /Tu- The initial wave packets <pi can be chosen optimally to minimize 
redundancy in the short-time correlations, and to obtain maximal information in a specific basis or for wave function 
structure in a given subspace of the original Hilbert space. 
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